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Abstract 

Oscillator models are central to the study of system properties such as entrainment or 
synchronization. Due to their nonlinear nature, few system-theoretic tools exist to analyze 
those models. The paper develops a sensitivity analysis for phase response curves, a funda- 
mental one-dimensional phase reduction of oscillator models. The proposed theoretical and 
numerical robustness analysis tools are illustrated on several models arising in the systems 
biology of cellular rhythms. 

1 Introduction 

Oscillator models — whose steady-state behavior is periodic rather than constant — are funda- 
mental to rhythmic modeling and they appear in many areas of engineering, physics, chemistry, 
and biology [THS]. 

Many oscillators are, by nature, open dynamical systems, that is, they interact with their 
environment [7]. As clocks or rhythm generators, their function often lies in their robust ability 
to respond to a particular input (entrainment) or to behave collectively in a network (synchro- 
nization or clustering). 

Starting with the pioneering work of Winfree [HE], the Phase Response Curve (PRC) of an 
oscillator has emerged as a fundamental input-output characteristic of oscillators. Analogously 
to the static gain of a transfer function, the PRC measures the steady-state (asymptotic) phase- 
shift caused by an impulsive input. Because it is phase dependent, the PRC is a curve rather 
than a scalar. In many situations, the PRC can be determined experimentally and provides 
unique data for the model identification of the oscillator. Likewise, numerical methods exist to 
compute the PRC from a state-space model of the oscillator. Finally, the PRC is the funda- 
mental mathematical information required to reduce a n-dimensional state-space model to the 
one-dimensional (phase) center manifold of a hyperbolic limit cycle. 

Motivated by the prevalence of this input-output characteristic and by the growing interest 
for system-theoretic questions related to oscillators, the present paper extends standard con- 
cepts of system theory to the space of PRCs. Comparing systems with a proper metric has been 
central to systems theory (see e.g. Zames [9|fTU] . Vinnicombe [IT] . Georgiou [12] for exemplative 
milestones), offering novel frameworks for system identification and robustness analysis. Follow- 
ing those steps, we aim at equipping the space of PRCs with the right metrics (accounting for 
natural invariance properties) and with the right sensitivity analysis tools to provide numerical 
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and mathematical grounds to system identification and robustness analysis of oscillator models. 
Although classical in their definition, several of these tools appear to be novel, in particular in 
the context of biological applications. 

We illustrate the potential use of the proposed approach by reconsidering recent examples of 
oscillator analysis in systems biology and in neurodynamics, suggesting that the proposed tools 
are adequate to address concrete questions of biologists and offer novel insight into existing 
questions. A preliminary version of this work was presented in |13j . 

The paper is organized as follows. Section [2] reviews the notion of PRCs characterizing 
the input-output behavior of an oscillator model in the neighborhood of a stable limit cycle. 
Section [3] defines several relevant metrics on (nonlinear) spaces of PRCs induced by natural 
equivalence properties. Section [3] develops the sensitivity analysis for oscillators in terms of the 
sensitivity of its periodic orbit and of its PRC. Section [5] illustrates how those tools permit to 
address system-theoretic questions arising for biological systems: robustness analysis, system 
identification, and model classification. Appendix [A] provides the numerical tools required to 
turn the abstract developments into concrete algorithms. 



2 Input (Infinitesimal) Phase Response Curves as Input— Output 
Characteristics of Oscillators 

In this section, we motivate and review the use of PRCs as fundamental input-output char- 
acteristics of oscillators. We first recall basic definitions about periodic orbits in state-space 
models of open dynamical systems (see |14H17j for details). We then summarize the standard 
phase reduction procedure which concentrates the phase behavior information of n-dimensional 
state-space models into one-dimensional phase models characterized by its angular frequency, 
its PRC, and a measurement map (see pl] [T8H2"0] for details). 

2.1 State-Space Model: Periodic Orbit and (Asymptotic) Phase 

We consider open dynamical systems described by nonlinear (single-input and single-outpu10) 
time-invariant state-space models 

x = f(x) + g(x)u (la) 
y = h{x) (lb) 

where states x(t) evolve on some subset X C M n , and input and output values u(t) and y(t) 
belong to subsets WCR and ^CK, respectively. The vector fields f : X — >• R n and g : X — > R n , 
and the measurement map h : X — )• y support all the usual smoothness conditions that are 
necessary for existence and uniqueness of solutions. We write cp(t,xo,u) for the solution of the 
initial value problem x = f(x) + g(x)u with x(0) = xq. 

We assume that the zero-input system x = f(x) admits a locally hyperbolic stable periodic 
orbit 7 C X with period T (and corresponding angular frequency uj = 2ir/T). Picking an initial 
condition x^ on the periodic orbit 7, the entire invariant set is described by the locus of the 
(nonconstant) T-periodic solution (j>(t, Xq,0) = x 7 (t), that is, 

j := {x £ X : x = <fi(t, Xq, 0),t E [0, T)} , (2) 

where the period T > is the smallest positive constant such that (f)(t,XQ,0) = <p(t + T, Xq,0) 
for all t G R. 



1 For presentation convenience, we consider single-input and single-output systems. All developments are easily 
generalizable to multiple-input and multiple-output systems. 
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Figure 1: The asymptotic phase map : B(-y) — > S 1 associates to each point x q in the basin B(j) 
a scalar phase Q{x q ) = 6 on the unit circle S 1 such that lim i _ i . +00 \\4>(t, x q , 0) — (p(t, x p , 0)|| 2 = 
with x p = x^iO/ui). 



The 6asin 0/ attraction of 7 (or oscillator stable set) is the maximal open set from which 
the periodic orbit 7 attracts 

B(-y) := {x £ X : lim dist(<£(t, x , 0), 7) = 0} (3) 

t— >+oo 

where dist(x,7) := inf y <= 7 ||x — 7/ 1 1 2 is the distance from the point x € X to the set 7 C ^ based 
on the Euclidean norm || • H2 in W 1 . 

Because a periodic orbit 7 is a one-dimensional (center) manifold in X , it is homeomorphic 
to the unit circle S 1 and then naturally parametrized in terms of a single scalar phase. Any point 
x p £ 7 can be characterized by a scalar phase d £ S 1 that uniquely determines the position of 
the point x p along the periodic orbit 7. The smooth bijective phase map : 7 — >■ S 1 associates 
to each point x p on the periodic orbit its phase i9, such that, 

x p = x 7 (tf/w). (4) 

This mapping is constructed such that the image of the reference point Xq is equal to and the 
progression along the periodic orbit (in absence of perturbation) produces a constant increase 
in '&. 

For hyperbolic stable periodic orbit, the notion of phase can be extended to any point x q 
in the basin -6(7) by defining the concept of asymptotic phase. The asymptotic phase map 
: 0(7) — > S 1 associates to each point x q in the basin 6(7) its asymptotic phase 6, such that 

lira U(t,x q ,0) -<^,x^/u;),0)|| 2 = 0. (5) 

Again, this mapping is constructed such that the image of Xq is equal to and the progres- 
sion along any orbit in #(7) (in absence of perturbation) produces a constant increase in 6. 
This mapping reduces to the phase map for points on the periodic orbit. (Main notations are 
illustrated on Fig. [TJ) 

Level sets of the asymptotic phase map 0, that is, sets of all points in the basin of 7 with 
the same asymptotic phase, are termed isochrons. Formally, the isochron Xq associated to the 
asymptotic phase 9 is the set Xq := {x G -8(7) : 0(x) = 8}. Considering hyperbolic periodic 
orbits, isochrons are codimension-1 submanifolds (diffeomorphic to R n_1 ) crossing the periodic 
orbit transversally and foliating the entire basin of attraction [21 J. 

For presentation convenience, we introduce the map x 7 : S 1 — > 7 which associates to 
each phase 6 a point <P(6/uj,Xq,0) = x 7 (#) on the periodic orbit. This map corresponds to 
a reparametrization of the periodic solution x 7 (-). 
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The 27r-periodic steady-state solution x 7 and the angular frequency u) can be calculated by 
solving the Boundary Value Problem (BVP) 

(aPT(0) - -f{x~<{9)) = (6a) 
to 

x 7 (2vr) -x 7 (0) = (6b) 

99(x 7 (0),w) = (6c) 

(where the prime denotes the derivative with respect to 9). The boundary conditions are 
given by the periodicity condition (|6b|) which ensures the periodicity of the map x 7 (-) and the 
phase condition (|6c|) which anchors a reference position x 7 (0) along the periodic orbit (see |22j 
for details). Numerical algorithms to solve this BVP are reviewed in Appendix lAl 

2.2 Reduced Phase Model: Phase Response Curves 

In the weak perturbation limit, that is, for small input u(t) = eu(t) (with eCl), any solution 
(j)(t, xq,u) of the system which starts in the neighborhood of the hyperbolic stable periodic orbit 
7 stays in its neighborhood. The n-dimensional state-space model can then be approximated 
by a one-dimensional phase model (see [6| fT9]. [2Q] ) 

9 = u) + eq u (9)u (7a) 

y = He) (7b) 

where the phase variable 6(t) evolves on the unit circle S 1 . This phase model is fully character- 
ized by the triplet: angular frequency u > 0, response map q u : S 1 — > K, and measurement map 

h : s 1 -> y. 

The map q u '. — y 1R is called the (input) iTifinitesiTncil Phase Response Curve (iPRC). It 
associates to each phase 6 the directional derivative of the asymptotic phase map at the point 
:r 7 (#) on the periodic orbit 7 along the control vector field g(x" f (9)), that is, 

q u (0):=D@(x'(e))[g(x'(e))} (8) 

where 

DQ(x)[n] := hm . (9) 

t-»0 t 

This directional derivative can be computed as the usual inner product in M. n 

DQ{^{e))[g{x<{6))\ = (V x @(^(0)),g(^(9))) (10) 

where Va;0(x 7 (#)) is the gradient of the asymptotic phase map at the point x 7 (#). The map 
q x : S 1 — > M. n which associates to each phase 9 the gradient of the asymptotic phase map at 
the point x 7 (9), that is, 

q x (9):=V x e(^(0)), (11) 

is the (state) infinitesimal Phase Response Curve (iPRC). 

The (state) iPRC q x : S 1 -> R n can be calculated by solving the BVP 

q'M + -U^(9))*q x (9) = (12a) 

fe (2vr) - q x (0) = (12b) 
(q x (9),f(x^9)))-u J = (12c) 
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where (j!2bj) is the boundary condition imposing the periodicity of q x (') and (|12cj) is a normal- 
ization condition ensuring a linear increase at rate uj of the phase variable 9 along unperturbed 
trajectories. Numerical methods to solve this BVP as a by-product of the periodic orbit com- 
putation are presented in Appendix lAl 

The measurement map h : S 1 — > y associates to each phase 6 the output value h(6) = 
(h o x 7 )(#). Without lost of generality, in this paper, we assume that the measurement map h 
is equal to the asymptotic phase map for points in the basin 6(7). The map h then reduces 
to the identity function idgi : S 1 — > S 1 : 9 h-> 6. The phase model is then characterized only by 
oj and q u (-). 

The approximate representation {ui , q u (-) , h(-)} of oscillators shares similar characteristics 
with the static gain of transfer function representation of linear time-invariant systems. Both 
representations capture asymptotic properties of the impulse response. They are independent 
of the complexity of the internal state-space representation of the oscillators. Moreover, those 
informations are available experimentally. 

3 Metrics in the Space of Phase Response Curves 

To address system-theoretic questions in the space of PRCs, it is of interest to equip this space 
with the differential structure of a Riemannian manifold. The differential structure permits to 
define local sensitivities in the tangent space. The Riemannian structure is convenient to recast 
analysis problems in an optimization framework, providing for instance a notion of steepest 
descent. The Riemannian structure also provides a norm in the tangent space and a (geodesic) 
distance between PRCs. 

Because PRCs are signals defined on the unit circle, the most obvious Riemannian structure 
is provided by the infinite-dimensional Hilbert space of square-integrable signals 

n°:={q:q(-)eC 2 (S\R)} (13) 

equipped with the standard inner product 

(£(•). C0>:= / mW)dO (14) 
Js 1 

and the associated norm 

U(-)h ■= X^WlW- (15) 

For technical reasons which will be clarify latter, we additionally assume that the first 
derivative of considered signals is also square-integrable. We then restrict the signal space to 
Hilbert space 

n 1 := {q : q(-) € £ 2 (S\ R), </(•) € £ 2 (§\M)} . (16) 

The space H 1 is a linear subspace of T~L° and it inherits its inner product (|14p and its norm (|15|) . 
This linear space structure is convenient for calculations but it fails to capture important equiva- 
lence properties between PRCs. In many applications, it is not meaningful to distinguish among 
PRCs that are related by a scaling factor and/or a phase-shift. 

Scaling equivalence The actual magnitude of the input signal acting on the system is not 
always well-known. This uncertainty about the input magnitude induces an (inversely pro- 
portional) uncertainty about the phase response magnitude. In those cases, we require that 

qi ~ q2 3a > : q 2 = qia (17) 
where a denotes any scaling factor (and ~ symbolizes the equivalence property). 
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q </■ qa. 


q ~ qa 


q q o a 


A: H 1 


B: ^VM >0 




C: ^ 1 /Shift(S 1 ) 


D: ^ 1 /(Shift(S 1 ) x M >0 ) 



Table 1: Combining equivalence properties defines different quotient spaces for phase response 
curves. 

Phase-shifting equivalence The choice of a reference position (associated to the initial 
phase) along the periodic orbit is often artificial. In those cases, we require that 

qi ~ q 2 3d E Shift(S 1 ) : q 2 = qi o a (18) 

where a denotes any phase-shift operator. A phase-shift operator corresponding to a phase-shift 
value A € S 1 is defined as the diffeomorphism a : S 1 — > S 1 : 9 i— > (9 + A) mod 27r. 

Those equivalence relations lead to the abstract — yet useful — concept of quotient space. 
Each point of a quotient space is defined as an equivalence class of signals. Since these equiv- 
alence classes are abstract objects, they cannot be explicitly used in numerical computations. 
Algorithms on quotient space work instead with representatives of these equivalence classes. 

Combining or not equivalence properties (|17p and (|18[) . we end up with four infinite- 
dimensional spaces: one Hilbert space and three quotient spaces, respectively, denoted by Qx 
with X = A, B, C, or D (Tab. [1]). In the next four subsections, we will equip each space with 
an appropriate Riemannian metric. 

3.1 Metric on Hilbert space H 1 

The simplest space structure is Hilbert space Qa '■=% ■ Each point q in Qa represents a signal. 
The Riemannian metric on Qa reduces to the inner product 

g q ^ q ,C q )-=(^,C q ) (19) 

and the norm in the tangent space T q Q\ at q is then given by 

U q \U ■= y/gq&itq) = yJ(Z q ,£ q ) = Ugh- (20) 

Because the space Qa is a linear space structure, the shortest path between two elements 
qi and q 2 on Qa is the straight line joining those elements. The natural (geodesic) distance 
between two points q\ and q 2 on Qa is then given by 

dist(gi, q 2 ) ■= \\qi - q 2 \\ 2 - (21) 

3.2 Metric on the quotient space H 1 /M >0 

The space capturing the scaling equivalence (fT7|) is the quotient space Qb := 'H 1 /M>o- Each 
element q in Qb represents an equivalence class 

q = [q] ■- {qa:a> 0}. (22) 
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Those equivalence classes are rays (starting at 0) in the total space Qb := H 1 . 
The Riemannian metric on Qb 

9q(tq,(q) : =9q(tiq,(q), (23) 

is induced by the metric on Qb 

^,C,):=%#, (24) 



which is invariant by scaling along the equivalence classes. As a consequence, it induces a metric 
g q (-, ■) on Qb- The norm in the tangent space T 9 Qb at q is given by 

lie, 



Uq\\ q -= V^^) = iSn ■ ( 25 ) 

v IMI2 

A signal representation of the tangent space at q G Qb relies on the decomposition of T^Qb 
into its vertical and horizontal subspaces. The vertical subspace Vq is the subspace of T^-Qb 
that is tangent to the equivalence class [q] 

Vq = {q(3 : P G R}. (26) 

The horizontal space Hq must be chosen such that T^-Qb = Hq © Vq. We choose Tig as the 
orthogonal complement of Vq for the metric ~g~q(-, •), that is, 

Hq = {r ] e TqQ B : g^r,, q(3) = 0}. (27) 

The orthogonal projection P^"q of a vector r) G T^Qb onto the horizontal space Hq is given by 

(28) 

y gq{qfi,ql3) (q,q) 

The distance between two points gi and q<i on Qb is defined as 

dist( gl , g2 ):=cos- 1 f || _ ( ^.f_ 2) || ) (29) 
(see [23] for metric on the unit sphere). 
3.3 Metric on the quotient space H 1 /Shift(§ 1 ) 

The space capturing the phase-shifting equivalence (fT8|) is the quotient space Qc '■=% / Shift (S 1 ). 
Each element q in Qc represents an equivalence class 

q = [q] = {q o a : a G Shift(S 1 )}. (30) 

Those equivalence classes are closed one-dimensional curves (due to the periodicity of the shift) 
on the infinite-dimensional hypersphere of radius ||g||2 in the total space Qc := H 1 . 
The Riemannian metric on Qc 

9q(tqXq) : =9q(£q-,(q), (31) 

is induced by the metric on Qc 

5f(£?,C?):=£ ? ,C ? >, (32) 
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which is invariant by phase-shifting along the equivalence classes. As a consequence, it induces 
a metric g q (-,-) on Qrj. The norm in the tangent space T q Qc at q is given by 



U g \U-=yj9 q (t; q ,ti q ) = U< I h. (33) 

The vertical space V q is the subspace of TqQc that is tangent to the equivalence class [q] 

Vq = {q'P : /3 G M} (34) 

(where <f has to belong to /^(S^M) to ensure the regularity of Vq). We choose the horizontal 
space T-Lq as the orthogonal complement of Vq for the metric g q (-, •), that is, 

Hq = {r, € TgQc : fljfa,^) = °>- ( 35 ) 
The orthogonal projection P^r\ of a vector r\ G TqQc onto the horizontal space %q is given by 

P h n ._ _ _ 9q(jhm , _ _ Mia' (36) 

The distance between two points q\ and q% on Qc is defined as 

dist(g-i, q 2 ) := min \\q 1 - q 2 o cr |) 2 = ||g x - q 2 o cr* || 2 (37) 
o-eShift(s 1 ) 

where <r* denotes the phase-shift operator achieving this minimization. It corresponds to the 
phase-shift operator maximizing the circular cross-correlation 

cr* = argmax (qi,q 2 °4 (38) 
o-eShift(§ 1 ) 

3.4 Metric on the quotient space T^/O^x) x Shift (S 1 )) 

The space capturing both scaling and phase-shifting equivalences (fTT|) ~ (fl~8j) is the quotient space 
Qd : = %V(^>o x Shift(S 1 )). Each element q in Qd represents an equivalence class 

q=[q] = {(qoa)a : a > 0,a 6 Shift^ 1 )}. (39) 

Based on the individual geometric interpretation of both equivalence properties, those equiva- 
lence classes are infinite cones in the total space Qd := H l , that is, the union of rays that start 
at and go through the closed one-dimensional curve of phase-shifted signals. 
The Riemannian metric on Qd 

9q(Zq,£q) ■= VqitqXq) , (40) 

is induced by the metric on Qd 

^(£o,C,) = %#, (41) 



which is invariant by scaling and phase-shifting along the equivalence classes. As a consequence, 
it induces a metric g q (-,-) on Qd- The norm in the tangent space T^Qd at q is given by 

lie, 



U q \U ■■= v^fo.£*) = g-. (42) 
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The vertical space Vq is the subspace of Tg-Qn that is tangent to the equivalence class [q] 

^ = M+^2:A,A6»}- (43) 

It is the direct sum of vertical spaces for equivalence properties individually. We choose the 
horizontal space %q as the orthogonal complement of Vq for the metric g q (-,-), that is, 

Uq = {rj £ TqQ D : g^r), qfa + q'f3 2 ) = 0}. (44) 
The orthogonal projection P-Jfrj of a vector r\ G T^Qd onto the horizontal space Tiq is given by 

Pgr) -=V- - i-n —a x gft - - La -,a J fa ( 45 ) 

H 9q{<lPi,qPi) gq{qP2,qP2) 

= n-^4~ q -^. (46) 

(q,q) {q',q'} 

The distance between two points (?i and q 2 on Qd is defined as 

dist( gi ,g 2 ):= min cos" 1 ( ^4l^r ) (47) 

= cos- 1 ( ffr: } ) (48) 

where cr* denotes the phase-shift operator achieving this minimization. It corresponds to the 
phase-shift operator maximizing the circular cross-correlation (see (|38j)). 

4 Sensitivity Analysis for Oscillators 

Sensitivity analysis for oscillators has been widely studied in terms of sensitivity analysis of 
periodic orbits |24H27j . In contrast, we found no reference on sensitivity analysis of PRCs in 
biology. Late in the preparation of the present manuscript, we discovered the recent work |28j 
which studies the sensitivity analysis of PRCs in the context of electronics. The developments 
in this section are closely related to those in [28J yet have been carried out independently. 

We summarize the sensitivity analysis for oscillators described by nonlinear time-invariant 
state-space models with one paramete: 
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x = f(x,X) + g(x,X)u (49a) 
y = h(x, X) (49b) 

where the constant parameter A belongs to some subset AC1. We first develop the sensitivity 
analysis of its periodic orbit and then provide the sensitivity analysis of its (i)PRC. 



4.1 Sensitivity Analysis of a Periodic Orbit 

The periodic orbit 7 of an oscillator model is characterized by its angular frequency u which 
measures the 'speed' of a solution along the orbit and by the 27r-periodic steady-state solu- 
tion xP{-) which describes the locus of this orbit in the state space. The sensitivity of both 
characteristics is important. 

2 For presentation convenience, we consider systems with a one-dimensional parameter space. All developments 
are easily generalizable to systems with a g-dimensional parameter space. 
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Given a nominal parameter value Ao , the angular frequency sensitivity is the scalar € 
defined as 



St. 



~dX 



Ao 



lim 

h->0 



h 



(50) 



where the notation *| A emphasizes the parameter value A at which the model characteristic * is 
evaluated. Likewise, the In -"periodic steady-state solution sensitivity is the 27r-periodic function 
Z x : S 1 -> R n defined as 



Z J.) ■- **L(.) 
W dX W 



r7M 



Ao 



lim • 



lA +/i 



IA 



where the explicit dependence of the 2-7r-periodic steady-state solution in A is given by 

* 7 (-)|a = 0(-/Ha,^q|a,O,A). 
From ([6|), we have, taking derivatives with respect to A, 

4(0) - -A(6)Z X (9) + 4^(0)^ - -b(9) = 

Z x {2n) - Z x (0) = 

<#rZaj(0) + <p w S u + (f\ = 

where we use the following short notations 



MP) 

b(9) 


:=|(^),A), 


fx 


■ 9 V ! 7 \\ 


v(6) 


:=f(^(e),x), 


fx 


7 \\ 

dX^o^^)- 



(51) 
(52) 

(53a) 

(53b) 
(53c) 

(54) 
(55) 
(56) 



4.2 Sensitivity Analysis of a Phase Response Curve 

The (input) infinitesimal phase response curve sensitivity is the 27r-periodic function Z qu : S 1 
R defined as 



Ao 



lim 



9«(0lAo+ft - ?«(' 



lAo 



From ([8]), we have, taking derivatives with respect to A, 

9g , 



dg , 



z qu {-) = (z qx (-),g(xH-),x)) + &(■), ^(-),x)z x (.) + ^0F(-),A) 



<9x 



<9A 



(57) 



(58) 



where the 27r-periodic function Z qx : S 1 — > M. n is the (state) infinitesimal phase response curve 
sensitivity defined as 



Z, 



<]x 1 



dq x 
dX 



Ao 



lim • 



9x(0L+ft- 



Ao 



From (|12p . we have, taking derivatives with respect to A, 



z' qx (6) + -A(eyz qx (e) + -c{ey qx {9) = o 

CO iu 

Z qx (27T)-Z qx (0) = 

(z qx (e),v(e)) + (q x (e), Zv(6)) -s^ = o 



(59) 



(60a) 

(60b) 
(60c) 
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where we use the following short notation 



<* W == g a^(^)> + A ) " hw,™ x)s - (61) 

Z,(9) :=^(^(e),X)Z i (9) + ^(Sr'(e),X). (62) 
4.3 Numerics of Sensitivity Analysis 

Numerical algorithms to solve BVPs (|53p and (|60p are reviewed in Appendix We stress that 
existing algorithms that compute periodic orbits and iPRCs are easily adapted to compute their 
sensitivity curves, essentially at the same numerical cost. 



5 Illustrations for Biological Systems 

In this section, we illustrate our approach on three system-theoretic questions arising for biolog- 
ical systems, emphasizing the novel insight provided by the proposed approach with respect to 
the existing literature. In the first illustration, we analyze the parametric robustness of a circa- 
dian oscillator model based on the sensitivity of its PRC. In the second illustration, we identify 
the parameter values of a simple circadian oscillator model in order to fit its experimental PRC. 
In the third illustration, we classify neural oscillator models based on their PRC. All numerical 
tests have been obtained with a MATLAB numerical code available from the authors. 



5.1 Parametric Robustness Analysis: a Case Study in a Quantitative Circa- 
dian Oscillator Model 

Testing the robustness of a model against parameter variations is a basic system-theoretic 
question. In a number of situations, the very purpose of modeling is to identify those parameters 
that influence a given system property. In the case of oscillator models, we propose scalar 
robustness measures to quantify the sensitivity of the angular frequency and of the iPRC to 
parameters. We apply those measures to a model of the circadian rhythm. 



5.1.1 Scalar Robustness Measure in the Space of PRCs 

The angular frequency u is a positive scalar number. The sensitivity of u> with respect to the 
parameter A is then also a scalar number S u , leading to the robustness measure defined as 

■= \S U \ (63) 

where | • | denotes the real absolute value function. 

In contrast, the iPRC (or its equivalence class) q belongs to a (nonlinear) space Q. The 
sensitivity of q is then a vector S q which belongs to the tangent space T q Q at q. A scalar 
robustness measure R q is defined as 

R q := \\S q \\ q = ^Jg q (S q ,S q ) (64) 

where ||-|| denotes the norm induced by the Riemannian metric g q (■, •) at q. It is the natural 
extension of robustness measures to (nonlinear) spaces. 
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When Q is a quotient space, the element q and the tangent vector S q are abstract objects. 
The evaluation of the robustness measure relies on the sensitivity Sq of the signal q defining the 
equivalence class in the total space 



R q 



pho_ 



^(pf^PfS^ (65) 



where P£ is the projection operator onto the horizontal space. The projection removes the 
component of the sensitivity which is tangent to the equivalence class. 

When analyzing a model with several parameters (A G A C M q ), all robustness measures 
R x (where x stands for any characteristic of the oscillator) collect in a g-dimensional vector the 
scalar robustness measure corresponding to each parameter. This vector is often normalized as 
follows 

R* = TT^T (66) 
\\Rx\\oo 

where || • ||oo denotes the maximum norm such that elements of R x belongs to the unit inter- 
val [0, 1]. This measure allows to rank model parameters according to their ability to influence 
the characteristic x. 

In addition, it is often meaningful to compute relative sensitivity measures, that is, relative 
changes in the model characteristic to relative changes in parameter values. 



5.1.2 Quantitative Circadian Oscillator Model 

In the literature, robustness analysis of circadian rhythms mostly studies the zero-input steady- 
state behavior (period, amplitude of oscillations, etc.) [29T31J and (empirical) phase-based per- 
formance measures |32H34| . 

Here, we illustrate our parametric robustness analysis on a quantitative circadian rhythm 
model. The Leloup-Goldbeter model is a computational model for the mammalian circadian 
clock based on the interlinked positive and negative feedback loops [55] . It describes the regu- 
latory interactions between the products of several genes (Per, Cry, and Bmall). This model 
possesses 16 state variables and 52 parameters. State-space model equations and nominal pa- 
rameter values are available in [351 Supporting Text]. It has been extensively studied through 
unidimensional bifurcation analyses and various numerical simulations of entrainment [35|.l36j . 

We develop our robustness analysis in the space Qd incorporating both scaling and phase- 
shifting equivalence properties. This is motivated by the uncertainty about the exact magnitude 
of the light input on the circadian oscillator and by the absence of precise experimental state 
trajectories preventing from defining a precise reference position (corresponding to the initial 
phase). 



5.1.3 Results 

A two-dimensional (R^, Rq u ) scatter plot in which each point corresponds to a model parameter 
reveals the shape and strength of the relationship between both angular frequency and iPRC 
robustness measures (Fig. [2]). In order to enlighten architecture rather than single-parameter 
properties, we grouped model parameters according to the mRNA path to which they belong 
(Per, Cry, or Bmall paths). In addition, we gathered parameters associated to interlocked 
paths in a last group. 

Firstly, the scatter plot reveals that most parameters exhibit both low angular frequency and 
iPRC sensitivities (points close to the origin); only few parameters display either high angular 
frequency or iPRC sensitivities. 
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Figure 2: Normalized robustness measures of the angular frequency R w and the iPRC R qu reveal 
the distinct sensitivity of three distinct gene paths ( Cry, Per, and Bmall ) . The three dashed 
lines are regression over the parameters of the three gene paths. Cry-path parameters exhibit 
low angular frequency and high iPRC sensitivities. Three of those parameters (v m c, v s c, and 
Kac) are more sensitive than others. 

Secondly, at a finer level of analysis, the scatter plot enables the identification of three 
distinct tendencies, one for each group of parameters corresponding to a particular mRNA 
path. Each tendency is materialized by a least-square linear regression on points belonging to 
a parameter group: 

• Cry-path: low frequency and high iPRC sensitivities; 

• Per-path: lower frequency than iPRC sensitivities; 

• Bmall -path: higher frequency than iPRC sensitivities. 

Moreover, among each tendency, some parameters are more sensitive than other. In par- 
ticular, three parameters associated to the Cry-path exhibit (very) high iPRC sensitivities and 
low angular frequency sensitivities: v m c, v s c, and Kac- Two of those parameters have been 
identified by numerical simulations as important for entrainment properties of the model with- 
out affecting the period (Kac in [35J and v m c in |36j). Our approach supports the importance 
of those two parameters and identifies the potential importance of a third one (v s c)- 

We stress that the conclusions in [35 , 36] rely on extensive simulations of the model to sim- 
ulate entrainment conditions while varying one parameter at a time. In contrast, the proposed 
analysis is systematic and computationally cheap. The plot in Fig. [2] is generated in less than 
a minute with our MATLAB code. 
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5.2 Parametric System Identification: a Case Study in a Qualitative Circa- 
dian Oscillator Model 

System identification deals with the problem of building mathematical models of dynamical 
systems based on observed data from the systems. In particular, parametric system identifi- 
cation aims at finding a set of parameter values in agreement with observed data for a given 
state-space model structure. 

In the literature, parameter values for circadian rhythm models are often determined em- 
pirically by trial-and-error methods due to the few pieces of experimental information about 
parameter values. 

In this section, we aim at providing a gradient-descent algorithm to identify a set of param- 
eter values which gives a PRC close to an experimental PRC (in our metric). We illustrate this 
algorithm on a qualitative circadian oscillator model. 



5.2.1 Gradient-descent algorithm in the space of PRCs 

A standard method to tackle the system identification problem is to recast it into an optimization 
framework. The minimization of an empirical cost V(A) yields the parameter estimate 

A = arg min y(A) (67) 
AeA 

where V(A) : A — > M>o penalizes the discrepancy between observed data from the system and 
prediction from the model. Local minimization is usually achieved with a gradient-descent 
algorithm requiring the computation of the gradient Va^(A). 

Given an 'experimental' PRC (or equivalence class) qo, a natural cost function V(X) is 
defined as 

V(X):=V(q\ x ) = ^di S t(q\ x ,q ) 2 (68) 

where dist(-, •) is the distance in the (nonlinear) space Q. The gradient (in the parameter space 
A) of this cost function with respect to the parameter Xj is given by 

V Xj V(X) = g q (grad, V(q\\), (69) 

where grad q V(q\ x ) and ^ are elements in the tangent space T q Q. 

When Q is a quotient space, the evaluation of the gradient Va^V^A) relies on representatives 
in the total space 

V Xj V(X) =g w (grad q V(q\x),P£-¥f) ( 70 ) 



where V(q) = V([q]) for all q € 



5.2.2 Qualitative Circadian Oscillator Model 

We illustrate the system identification on a qualitative circadian rhythm model. The Goodwin 
oscillator is a cyclic feedback system where metabolites repress the enzymes which are essential 
for their own synthesis by inhibiting the transcription of the molecule DNA to messenger RNA 
(mRNA) [37j. It can be described as the cyclic interconnection of three first-order subsystems 
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and a monotone static nonlinearity 



Hi 




i + (ui/ey 



1 



for % = 1, . . . , 3 



where the cyclic interconnection is given by = — g/4, U2 = Ui, U3 = y2, and U4 = ii ex t + 2/3- A 
dimensionless form of this system is equivalent to impose K2 = K3 = t\ = 9 = 1. 

To simplify the analysis (but without loss of generality), we reduce the parameter space to 
two dimensions: we impose equal time-constants in H2 and #3 (T2 = T3 = r) and fix the Hill 
coefficient p = 20. This high coefficient is justified by the necessity to get periodic orbits (p > 8) 
and strong enough differences between iPRC shapes in the parameter space. The results for 
weaker coefficient are similar but less marked. The parameter space reduces to {K\,t) € K>q- 

An 'experimental-like' iPRC is chosen as the iPRC computed for a quantitative circadian 
rhythm model of Drosophila which has proven to agree with experimental data j38|.l39| . 

In this context, it is meaningful to perform the identification in the space <2d incorporating 
both scaling and phase-shifting equivalence properties for the same reasons as in the previous 
illustration. 

5.2.3 Results 

The Goodwin model exhibits stable oscillations in a C-shaped region of the reduced parameter 
space (Fig. 13(a)). The border of this region corresponds to a supercritical Andronov-Hopf 
bifurcation through which the model single equilibrium looses its stability. The contour levels 
of the cost function — which have been computed in the whole C-shaped region to make results 
interpretation easier — reveal two local minima and then the nonconvexity of the cost function. 

Picking initial guess values for model parameters, the gradient-descent algorithm minimizes 
the cost function following a particular path in the parameter space (Fig. [3(a)). The cost 
function value decreases at each step of the algorithm along this path (Fig. [3(b)). The optimal 
iPRC fits very well the 'experimental' one, in contrast to the initial iPRC (Fig. [3(c)). 

Due to the nonconvexity of the cost function, two paths starting from close initial points 
may evolve towards different local minima (Fig. [3(a)). However, the cost function is (almost) 
symmetric with respect to a unitary time-constant r and both local minima correspond to 
similar iPRCs (up to a scaling factor and a phase-shift). 

5.3 Parametric Model Classification: a Case Study on a Neural Oscillator 



Model classification groups model in classes which share common qualitative and/or quantitative 
characteristics. 

In the literature, models of neurons are often classified into two types based on the bifur- 
cation that gives birth to periodic firing |4L)| . Class-I excitable neurons arise from saddle-node 
on invariant circle bifurcations and can theoretically fire at arbitrarily low finite frequencies. 
Class-II excitable neurons arise from a subcritical or supercritical Andronov-Hopf bifurcations 
and possess a nonzero minimum frequency of firing. 

In this section, we compare this model classification to a classification based on the distance 
to canonical iPRCs in the space of iPRCs. 



Model 
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Figure 3: (a) The cost function (gray levels) between an experimental iPRC and the (input) 
iPRCs exhibits a nonconvex behavior in the reduced parameter space. The gradient-descent 
algorithm follows the path indicated by dots, (b) The distance along the path followed by the 
gradient-descent algorithm decreases with the iteration number, (c) The shape of the optimal 
iPRC is closer to the reference iPRC than the initial iPRC. 



5.3.1 Model Classification Scheme in the Space of PRCs 

A strong relationship between the bifurcation type and the shape of the iPRC has been demon- 
strated [20 , 401EET] . On the one hand, the iPRC of class-I excitable neurons near the bifurcation 
is nonnegative or nonpositive and approximated by 

9Sn(0) = — [1-cos(0-0 S n)] (71) 

where csn and #sn are model dependent constants. One the other hand, the iPRC of class-II 
excitable neurons near the bifurcation has both positive and negative parts and is approximated 
by 

q H (9) = CH sin(0 - On) (72) 

- U>H\ 

where ch, ^h, and wh are model dependent constants. 

We propose to classify models in the parameter space based on the distance between the 
model iPRC and canonical iPRCs 

^ ^ I class-I if dist(g,gsN) < dist(g,g H ) ^ 
1 class-II if dist(g, qsn) > dist(g, gn) 

where dist(-, •) is the distance in the space Q. 
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5.3.2 Neuron Oscillator Model 

The Morris-Lecar model is a popular two-dimensional reduced model of excitable neurons [H] 



CV 



w = 



-5Ca^oo(^)(F - V Ca ) " 9^w{V ~ V K ) 

gdv- v h )+i app 

Wqo(V) - w 

r w {V) 



(74) 
(75) 



(76) 



where 



"loo 00 
™oo00 



0.5(1 +ianh((V -Vi)/V 2 )) 
0.5(1 +tanh((V-F 3 )/V 4 )) 



(77) 
(78) 
(79) 



and 



r w {V) 



l/cosh((y-F 3 )/(2T4)). 



(80) 



The applied current I app plays the role of input. 

This model exhibits both classes of excitability for different parameter values [43, 44J . In 
particular, the calcium conductance ~g~Q a modifies the nature of the bifurcation giving birth to 
the periodic orbit. For large values of g Ca , the model exhibits a class-I excitability (saddle- 
node on invariant circle bifurcation). For smaller value of g Ca , the model exhibits a class-II 
excitability (Andronov-Hopf bifurcation). 

In this context, it is meaningful to classify model based on a distance in the space Qd 
incorporating both scaling and phase-shifting equivalence properties. We are indeed interested 
only in comparing the qualitative shape of iPRCs. 

5.3.3 Results 

Standard classification scheme based on the excitability type defines an horizontal separation 
in the two-dimensional parameter space (-f app , ffca) (Fig- SKa))- 

The classification scheme based on the iPRC shape gives a very different separation in the 
parameter space (Fig. HJb)). This classification scheme allows one neuron (for one value of g Ca ) 
to pass from one class to another (crossing of the separation) for different values of applied 
current I app . iPRCs computed for several points close to the bifurcation boundary confirm the 
classification based on the qualitative shape of iPRCs. 

For class-II oscillators, we observe that the correspondence between the bifurcation-based 
classification and the PRC-based classification is limited to a narrow region in the neighborhood 
of the bifurcation. 

6 Conclusion 

The paper proposes sensitivity tools to analyze oscillators as open dynamical systems. We 
showed that, under the weak perturbation assumption, state-space models can be reduced to 
phase models characterized by their angular frequencies and their (input infinitesimal) phase 
response curves. 

We proposed to base metrics in the space of dynamical systems on metrics in the space 
of (input infinitesimal) phase response curves. We focussed on nonlinear spaces induced by 
equivalence properties among phase response curves. 
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Figure 4: Bifurcation-based and PRC-based classifications determine very different subsets of 
the parameter space for the Morris-Lecar model. Parameter values: C = 20 /xF/cm 2 , g K = 
8 mS/cm 2 , g L = 2 mS/cm 2 , F Ca = 120 mV, V K = -80 mV, V L = -60 mV, Vi = -1.2 mV, 
V 2 = 18 mV, y 3 = 12 mV, V 4 = 17.4 mV, <j) = 1/15 s _1 . 

We also introduced the sensitivity analysis of periodic orbits and their phase response curves. 

The combination of the metrics and the sensitivity analysis allow to answer system-theoretic 
questions. We illustrated our approach on system-theoretic questions arising for biological sys- 
tems. While preliminary, those illustrations suggest that the proposed approach is numerically 
efficient and may provide novel insight in several questions of interest for oscillator modeling. 

A Numerical Tools 

Several numerical algorithms exist for the numerical computation of periodic orbits 22, lol[~R)1. 
Most algorithms recast the periodic orbit computation as a two-point boundary value problems. 
Numerical boundary value methods fall into two classes: 

1. shooting methods generate trajectory segments using a numerical time integration and 
match segment end points with each other and the boundary conditions; 

2. global methods project the differential equations onto a finite dimensional space of discrete 
closed curves that satisfy the boundary conditions. 

Both methods yield a set of (nonlinear) equations that are solved with root finding algorithms, 
usually Newton's method. 

In this appendix, we summarize popular algorithms for the computation of periodic orbits. 
Then we emphasize how the computation of the (state) infinitesimal phase response curve is a 
cheap by-product of this computation. Finally we extend those algorithms for the computation 
of oscillator sensitivities: angular frequency, steady-state periodic solution, and infinitesimal 
phase response curve sensitivities. 



18 



Gi 
Hi 

Gi 
Hi 

P 

Ef 



Forward Shooting 



Trapezoidal Scheme 



4> a) -x i+ i 

In. 



hi d± I hi_ ~ n \ 



d± (hi f. n \Y 
1 r 

]V+T J (AT+l)n 

(hi r- o \ 



Xl 2 w 



[/ (sc i ,A) + /(5 i+ i,A)] 



-/n + ^/,(x m ,A) 

-1^ [/(£«, A) + A)] 



l^f/a: A)* 



J_ Hino- I h 0+ h l fejv-i+fejy fejy 

2vr Ulcl & 2' 2 2 '2 



: % + C {9i+i)* qi+i] 



Table 2: Residual maps ri(xs,ui), linear block elements Gi, Hi, and hi, and adjoint linear block 
elements Gi and Hi for two one-step numerical algorithms (i = 0, 1, . . . ,N — 1). 



A.l Numerical Computation of Periodic Orbits 

A periodic orbit 7 is characterized by the 2-7r-periodic steady-state solution x 7 : 8 1 — > 7 describ- 
ing a closed curve in the state space and the angular frequency u € M>o (or equivalently the 
period T) which solve the boundary value problem ([6]). 

Introducing a (nonuniform) partition S of the unit circle S 1 

5 : = 6 < 61 < ■ ■ ■ < 6 N = 2vr, (81) 

the 2-7r-periodic steady-sate solution x 1 is numerically approximated by a closed discrete curve 
in the state space X. A discrete curve is a set of points {xq,x\, . . . ,xn} associated to the 
set of phases such that Xi approximates x 7 (#j) for all i = 0, 1, ... ,N. This discrete curve is 
closed, that is, x n = xq, which reflects the periodicity of the solution x 7 . In the following, 
the circle partition 5 is fixed and the discrete curve is numerically represented by the vector 
x$ = [xq ,x]- , . . . ,xJj] T . We denote hi = 9i + \ — Q, L . 

Equations for approximate periodic orbits take then the form of N n-dimensional vector 
equations 

n{xs,u) = 0, i = 0,l,...,N -1, (82) 

where different residual maps lead to different numerical methods (Tab.[2])- Those equations 
are completed by the periodicity condition 

r N (x s ,uj) := x N - x = (83) 

and the phase condition 

r v (xs,u) : = <p(xs,u}, A) = (84) 

where (p is a discretized version of the phase condition (p. 

We solve this set of (nonlinear) equations r(xg,u)) = with the root finding Newton's 

method. Starting from an initial guess fxjj°\ c*/°)J , this method iteratively update the solution 
xf +1) =xf ) +Axf ) and J fe+1 ) = u,« + Aw< fc >. (85) 



19 



Correction terms are computed by solving the linear problem 
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where A has a particular block structure for one-step schemes and b, c, and d are also defined 
by blocks 
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Expressions of block elements Gi, Hi, and bi depends on the methods used to generate residual 
maps ri(xs, oj) = 0, with i = 0, 1, . . . , N — 1, for approximate periodic orbits (Tab. [2]). 

The main computational effort in one iteration arises from the evaluation of the (N + l)n x 
(N + l)n structured matrix A whose block elements are computed through fundamental solution 
time integrations or Jacobian matrix evaluations. 



A. 2 Numerical Computation of Phase Response Curves 

The (state) infinitesimal phase response curve q x : S 1 — > W 1 of a periodic orbit is the solution 

of the boundary value problem ([12]) . 

The (state) infinitesimal phase response curve is numerically approximated by a set of points 

{go, Qi, ■ ■ ■ , Qn} associated to the phases in the partition 5 such that = go- 
Following the same procedure as for approximate periodic orbits, equations for approximate 

infinitesimal phase response curves take the form of (N + l)n linear equations 



Aq 5 = 

where the matrix A has the same structure as the matrix A 

Go —Hq 

A 



(89) 
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Block elements of A can be constructed based on numerical computations done for the periodic 
orbit computation (Tab. [2]). 

The matrix A is by construction singular with a simple rank deficiency. This rank deficiency 
is overcome by adding a normalization condition for q$. Using (|12cp . we have 
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where vg is the approximate tangent vector to the periodic orbit and P is a ponderation matrix 
which depends on the method class. We seek a system of defining equations that is square 
and regular. A standard way to tackle this issue is to border the matrix A as follows (see |22[ 
Theorem 5.8] for details) 
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with d q 7^ 0, cj 



vTP, and b q £ range(^4) (for example b q 
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A. 3 Numerical Computation of Oscillator Sensitivities 



The angular frequency sensitivity S u £ W- Xq and the 2-7r-periodic steady-sate solution sensitivity 
Z x : S 1 — > M nx<? are the solutions of the linear boundary value problem (|53p . Equations for 
approximate periodic orbit sensitivities take the form of a system of linear equations 
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(93) 



where = — and Ef depends on the numerical method used (Tab. [2]). 

The infinitesimal phase response curve sensitivity Z q : S 1 — > W ixq is the solution of the 
linear boundary value problem (|60p . Equations for approximate infinitesimal phase response 
curve sensitivities take the form of a system of linear equations 



A b q 
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(94) 



where = S w — qj PZ V ^ and E\ depends on the numerical method used (Tab. [2]). 

In both cases, square matrices in left hand sides are identical to the matrices used for the 
computation of the periodic orbit and the infinitesimal phase response curve, respectively. The 
only additional computation effort arises from the evaluation of the right hand side. 
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